home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cagddist.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  10.8 KB  |  365 lines
  1. /******************************************************************************
  2. * CagdDist.c - Distance computations.                          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, May 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. #define SELF_INTER_REL_EPS 100
  13.  
  14. static CagdPtStruct
  15.     *GlblSrfDistPtsList = NULL;
  16.  
  17. static void CagdSrfDistFindPointsAux(CagdSrfStruct *Srf, CagdRType Epsilon,
  18.                              CagdBType SelfInter);
  19. static void SrfDistAddZeroPoint(CagdRType u, CagdRType v, CagdRType Epsilon);
  20.  
  21. /******************************************************************************
  22. * Given a curve and a point, find the closest point (if MinDist) or the far   *
  23. * location (if !MinDist) from the curve to the given point.              *
  24. *   Returned is the parameter value of the curve.                  *
  25. ******************************************************************************/
  26. CagdRType CagdDistCrvPoint(CagdCrvStruct *Crv, CagdPType Pt, CagdBType MinDist,
  27.                                  CagdRType Epsilon)
  28. {
  29.     CagdRType TMin, TMax, t,
  30.     ExtremeDist = MinDist ? INFINITY : -INFINITY;
  31.     CagdPtStruct *TPt,
  32.     *Pts = CagdLclDistCrvPoint(Crv, Pt, Epsilon);
  33.  
  34.     /* Add the two end point of the domain. */
  35.     CagdCrvDomain(Crv, &TMin, &TMax);
  36.     TPt = IritMalloc(sizeof(CagdPtStruct));
  37.     TPt -> Pt[0] = TMin;
  38.     TPt -> Pnext = Pts;
  39.     Pts = TPt;
  40.     TPt = IritMalloc(sizeof(CagdPtStruct));
  41.     TPt -> Pt[0] = TMax;
  42.     TPt -> Pnext = Pts;
  43.     Pts = TPt;
  44.  
  45.     for (TPt = Pts, t = TMin; TPt != NULL; TPt = TPt -> Pnext) {
  46.     int i;
  47.     CagdPType EPt;
  48.     CagdRType
  49.         Dist = 0.0,
  50.         *R = CagdCrvEval(Crv, TPt -> Pt[0]);
  51.  
  52.     CagdCoerceToE3(EPt, &R, - 1, Crv -> PType);
  53.  
  54.     for (i = 0; i < 3; i++)
  55.         Dist += SQR(EPt[i] - Pt[i]);
  56.  
  57.     if (MinDist) {
  58.         if (Dist < ExtremeDist) {
  59.         t = TPt -> Pt[0];
  60.         ExtremeDist = Dist;
  61.         }
  62.     }
  63.     else {
  64.         if (Dist > ExtremeDist) {
  65.         t = TPt -> Pt[0];
  66.         ExtremeDist = Dist;
  67.         }
  68.     }
  69.     }
  70.  
  71.     return t;
  72. }
  73.  
  74. /******************************************************************************
  75. * Given a curve and a point, find the local extremum distance points on the   *
  76. * curve to the given point. Return is a list of parametr value with local     *
  77. * extremum.                                      *
  78. *   Computes the zero set of (Crv(t) - Pt) . Crv'(t).                  *
  79. ******************************************************************************/
  80. CagdPtStruct *CagdLclDistCrvPoint(CagdCrvStruct *Crv, CagdPType Pt,
  81.                                  CagdRType Epsilon)
  82. {
  83.     int i;
  84.     CagdCrvStruct *ZCrv,
  85.     *DCrv = CagdCrvDerive(Crv);
  86.     CagdPType MinusPt;
  87.     CagdPtStruct *ZeroSet;
  88.  
  89.     for (i = 0; i < 3; i++)
  90.     MinusPt[i] = -Pt[i];
  91.  
  92.     Crv = CagdCrvCopy(Crv);
  93.     CagdCrvTransform(Crv, MinusPt, 1.0);
  94.  
  95.     ZCrv = CagdCrvDotProd(Crv, DCrv);
  96.     CagdCrvFree(Crv);
  97.     CagdCrvFree(DCrv);
  98.  
  99.     ZeroSet = CagdCrvZeroSet(ZCrv, 1,  Epsilon);
  100.     CagdCrvFree(ZCrv);
  101.  
  102.     return ZeroSet;
  103. }
  104.  
  105. /******************************************************************************
  106. * Given a curve and a line, find the closest point (if MinDist) or the far    *
  107. * location (if !MinDist) from the curve to the given line.              *
  108. *   Returned is the parameter value of the curve.                  *
  109. ******************************************************************************/
  110. CagdRType CagdDistCrvLine(CagdCrvStruct *Crv, CagdLType Line,
  111.                       CagdBType MinDist, CagdRType Epsilon)
  112. {
  113.     CagdRType TMin, TMax, t,
  114.     ExtremeDist = MinDist ? INFINITY : -INFINITY;
  115.     CagdPtStruct *TPt,
  116.     *Pts = CagdLclDistCrvLine(Crv, Line, Epsilon);
  117.  
  118.     /* Add the two end point of the domain. */
  119.     CagdCrvDomain(Crv, &TMin, &TMax);
  120.     TPt = IritMalloc(sizeof(CagdPtStruct));
  121.     TPt -> Pt[0] = TMin;
  122.     TPt -> Pnext = Pts;
  123.     Pts = TPt;
  124.     TPt = IritMalloc(sizeof(CagdPtStruct));
  125.     TPt -> Pt[0] = TMax;
  126.     TPt -> Pnext = Pts;
  127.     Pts = TPt;
  128.  
  129.     for (TPt = Pts, t = TMin; TPt != NULL; TPt = TPt -> Pnext) {
  130.     CagdPType EPt;
  131.     CagdRType
  132.         Dist = 0.0,
  133.         *R = CagdCrvEval(Crv, TPt -> Pt[0]);
  134.  
  135.     CagdCoerceToE2(EPt, &R, - 1, Crv -> PType);
  136.  
  137.     Dist = EPt[0] * Line[0] + EPt[1] * Line[1] + Line[2];
  138.     Dist = ABS(Dist);
  139.  
  140.     if (MinDist) {
  141.         if (Dist < ExtremeDist) {
  142.         t = TPt -> Pt[0];
  143.         ExtremeDist = Dist;
  144.         }
  145.     }
  146.     else {
  147.         if (Dist > ExtremeDist) {
  148.         t = TPt -> Pt[0];
  149.         ExtremeDist = Dist;
  150.         }
  151.     }
  152.     }
  153.  
  154.     return t;
  155. }
  156.  
  157. /******************************************************************************
  158. * Given a curve and a line, find the local extremum distance points on the    *
  159. * curve to the given line. Return is a list of parametr value with local      *
  160. * extremum.                                      *
  161. *   Let Crv be x(t), y(t). By substituting x(t) and y(t) into the line          *
  162. * equation, we derive the distance function. Its zero set, combined with the  *
  163. * zero set of its derivative provide the needed extrema.              *
  164. ******************************************************************************/
  165. CagdPtStruct *CagdLclDistCrvLine(CagdCrvStruct *Crv, CagdLType Line,
  166.                                  CagdRType Epsilon)
  167. {
  168.     CagdPType Pt;
  169.     CagdCrvStruct *TCrv1, *CrvX, *CrvY, *CrvZ, *CrvW, *DistCrv, *DerivDistCrv;
  170.     CagdPtStruct *ZeroSet1, *ZeroSet2, *TPt;
  171.  
  172.     CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
  173.     if (CrvZ)
  174.     CagdCrvFree(CrvZ);
  175.  
  176.     Pt[0] = Pt[1] = Pt[2] = 0.0;
  177.     CagdCrvTransform(CrvX, Pt, Line[0]);
  178.     CagdCrvTransform(CrvY, Pt, Line[1]);
  179.     TCrv1 = CagdCrvAdd(CrvX, CrvY);
  180.     CagdCrvFree(CrvX);
  181.     CagdCrvFree(CrvY);
  182.     if (CrvW) {
  183.     CagdCrvTransform(CrvW, Pt, Line[2]);
  184.     DistCrv = CagdCrvAdd(TCrv1, CrvW);
  185.     CagdCrvFree(CrvW);
  186.     CagdCrvFree(TCrv1);
  187.     }
  188.     else {
  189.     Pt[0] = Line[2];
  190.     CagdCrvTransform(TCrv1, Pt, 1.0);
  191.     DistCrv = TCrv1;
  192.     }
  193.  
  194.     ZeroSet1 = CagdCrvZeroSet(DistCrv, 1,  Epsilon);
  195.  
  196.     DerivDistCrv = CagdCrvDerive(DistCrv);
  197.     CagdCrvFree(DistCrv);
  198.  
  199.     ZeroSet2 = CagdCrvZeroSet(DerivDistCrv, 1,  Epsilon);
  200.     CagdCrvFree(DerivDistCrv);
  201.  
  202.     if (ZeroSet1 == NULL)
  203.     return ZeroSet2;
  204.     else if (ZeroSet2 == NULL)
  205.     return ZeroSet1;
  206.     else {
  207.     for (TPt = ZeroSet1; TPt -> Pnext != NULL; TPt = TPt -> Pnext);
  208.  
  209.     TPt -> Pnext = ZeroSet2;
  210.  
  211.     return ZeroSet1;
  212.     }
  213. }
  214.  
  215. /******************************************************************************
  216. * Given two curves, creates a scalar surface representing the distance square *
  217. * between them.                                      *
  218. ******************************************************************************/
  219. CagdSrfStruct *CagdSrfDistCrvCrv(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
  220. {
  221.     CagdSrfStruct *TSrf, *DiffSrf, *DotProdSrf,
  222.     *Srf1 = CagdPromoteCrvToSrf(Crv1, CAGD_CONST_U_DIR),
  223.     *Srf2 = CagdPromoteCrvToSrf(Crv2, CAGD_CONST_V_DIR);
  224.  
  225.     if (CAGD_IS_BSPLINE_SRF(Srf1) || CAGD_IS_BSPLINE_SRF(Srf2)) {
  226.     CagdRType UMin1, UMax1, VMin1, VMax1, UMin2, UMax2, VMin2, VMax2;
  227.  
  228.     if (CAGD_IS_BEZIER_SRF(Srf1)) {
  229.         TSrf = CnvrtBezier2BsplineSrf(Srf1);
  230.         CagdSrfFree(Srf1);
  231.         Srf1 = TSrf;
  232.     }
  233.     if (CAGD_IS_BEZIER_SRF(Srf2)) {
  234.         TSrf = CnvrtBezier2BsplineSrf(Srf2);
  235.         CagdSrfFree(Srf2);
  236.         Srf2 = TSrf;
  237.     }
  238.  
  239.     CagdSrfDomain(Srf1, &UMin1, &UMax1, &VMin1, &VMax1);
  240.     CagdSrfDomain(Srf2, &UMin2, &UMax2, &VMin2, &VMax2);
  241.     BspKnotAffineTrans(Srf1 -> VKnotVector,
  242.                Srf1 -> VLength + Srf1 -> VOrder,
  243.                VMin2 - VMin1, (VMax2 - VMin2) / (VMax1 - VMin1));
  244.     BspKnotAffineTrans(Srf2 -> UKnotVector,
  245.                Srf2 -> ULength + Srf1 -> VOrder,
  246.                UMin1 - UMin2, (UMax1 - UMin1) / (UMax2 - UMin2));
  247.     }
  248.  
  249.     DiffSrf = CagdSrfSub(Srf1, Srf2);
  250.     DotProdSrf = CagdSrfDotProd(DiffSrf, DiffSrf);
  251.     
  252.     CagdSrfFree(Srf1);
  253.     CagdSrfFree(Srf2);
  254.     CagdSrfFree(DiffSrf);
  255.  
  256.     return DotProdSrf;
  257. }
  258.  
  259. /******************************************************************************
  260. * Given a scalar surface representing the distance square between two curves, *
  261. * find the zero set of the surface, if any and return it.              *
  262. *   The given surface is a non negative surface and zero set is its minima.   *
  263. *   The returned points will contail the two parameter values of the two      *
  264. * curves that intersect in point's X & Y.                      *
  265. ******************************************************************************/
  266. CagdPtStruct *CagdSrfDistFindPoints(CagdSrfStruct *Srf, CagdRType Epsilon,
  267.                             CagdBType SelfInter)
  268. {
  269.     GlblSrfDistPtsList = NULL;
  270.  
  271.     if (Srf -> GType == CAGD_SBEZIER_TYPE) {
  272.     Srf = CnvrtBezier2BsplineSrf(Srf); /* To keep track on the domains. */
  273.     CagdSrfDistFindPointsAux(Srf, Epsilon, SelfInter);
  274.     CagdSrfFree(Srf);
  275.     }
  276.     else
  277.     CagdSrfDistFindPointsAux(Srf, Epsilon, SelfInter);
  278.  
  279.     return GlblSrfDistPtsList;
  280. }
  281.  
  282. /******************************************************************************
  283. * Aux. function of CagdSrfDistFindPoints - does the hard work.              *
  284. ******************************************************************************/
  285. static void CagdSrfDistFindPointsAux(CagdSrfStruct *Srf, CagdRType Epsilon,
  286.                              CagdBType SelfInter)
  287. {
  288.     int i, j,
  289.     ULength = Srf -> ULength,
  290.     VLength = Srf -> VLength;
  291.     CagdRType
  292.     *XPts = Srf -> Points[1];
  293.  
  294.     for (i = 0; i < ULength; i++) {
  295.     for (j = 0; j < VLength; j++) {
  296.         if (*XPts++ <= 0.0) {
  297.         CagdRType UMin, UMax, VMin, VMax;
  298.  
  299.         CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  300.  
  301.         if (SelfInter) {
  302.             CagdRType
  303.             SelfInterEps = Epsilon * SELF_INTER_REL_EPS;
  304.  
  305.             if (UMax - UMin < SelfInterEps &&
  306.             VMax - VMin < SelfInterEps &&
  307.             UMax - VMin < SelfInterEps &&
  308.             VMax - UMin < SelfInterEps)
  309.                 return;
  310.         }
  311.  
  312.         /* The surface may have a zero here. Test the domain size   */
  313.         /* to make sure it makes sense to subdivide.            */
  314.         if (UMax - UMin < Epsilon && VMax - VMin < Epsilon) {
  315.             SrfDistAddZeroPoint((UMin + UMax) / 2.0,
  316.                     (VMin + VMax) / 2.0, Epsilon);
  317.         }
  318.         else {
  319.             /* Subdivide the surface and invoke recursively. */
  320.             CagdSrfStruct *Srf1, *Srf2;
  321.             CagdRType t;
  322.             CagdSrfDirType Dir;
  323.  
  324.             if (UMax - UMin > VMax - VMin) {
  325.             t = (UMin + UMax) / 2.0;
  326.             Dir = CAGD_CONST_U_DIR;
  327.             }
  328.             else {
  329.             t = (VMin + VMax) / 2.0;
  330.             Dir = CAGD_CONST_V_DIR;
  331.             }
  332.             Srf1 = CagdSrfSubdivAtParam(Srf, t, Dir);
  333.             Srf2 = Srf1 -> Pnext;
  334.             CagdSrfDistFindPointsAux(Srf1, Epsilon, SelfInter);
  335.             CagdSrfDistFindPointsAux(Srf2, Epsilon, SelfInter);
  336.             CagdSrfFree(Srf1);
  337.             CagdSrfFree(Srf2);
  338.         }
  339.         return;
  340.         }
  341.     }
  342.     }
  343. }
  344.  
  345. /******************************************************************************
  346. * Aux. function of CagdSrfDistFindPoints - does the hard work.              *
  347. ******************************************************************************/
  348. static void SrfDistAddZeroPoint(CagdRType u, CagdRType v, CagdRType Epsilon)
  349. {
  350.     CagdPtStruct *Pt;
  351.  
  352.     for (Pt = GlblSrfDistPtsList; Pt != NULL; Pt = Pt -> Pnext) {
  353.     if (ABS(Pt -> Pt[0] - u) < Epsilon * 10 &&
  354.         ABS(Pt -> Pt[1] - v) < Epsilon * 10)
  355.         return;                     /* Point is already there. */
  356.     }
  357.  
  358.     Pt = IritMalloc(sizeof(CagdPtStruct));
  359.  
  360.     Pt -> Pt[0] = u;
  361.     Pt -> Pt[1] = v;
  362.     Pt -> Pnext = GlblSrfDistPtsList;
  363.     GlblSrfDistPtsList = Pt;
  364. }
  365.